function Panorama( src1,src2,RANSAC_distance,RANSAC_ratio )
%Step1 load images
img_left_origin=imread(src1);
img_right_origin=imread(src2);
if size(size(img_left_origin),2)==3
    img_left=single(rgb2gray(img_left_origin));
else 
    img_left=single(img_left_origin);
end
if size(size(img_right_origin),2)==3
    img_right=single(rgb2gray(img_right_origin));
else
    img_right=single(img_right_origin);
end
% imshow(img_left,[0,255]);
% imshow(img_right,[0,255]);

%Step2 detect feature points
peak_thresh=0.1;
edge_thresh=1.5;
[f_left,d_left]=vl_sift(img_left,'PeakThresh', peak_thresh,'edgethresh', edge_thresh);
[f_right,d_right]=vl_sift(img_right,'PeakThresh', peak_thresh,'edgethresh', edge_thresh);

%Step3 compute distances
size_left=size(d_left);
size_right=size(d_right);
distance_left_right=zeros(size_left(2),size_right(2));
for i=1:size_left(2)
    for j=1:size_right(2)
        minus=d_left(:,i)-d_right(:,j);
        distance_left_right(i,j)=sum(minus.*minus);
    end
end
distance_left_right_norm = 100*(distance_left_right - min(min(distance_left_right)))/(max(max(distance_left_right)) - min(min(distance_left_right)));

%Step4 select putatice matches
myNorm = distance_left_right_norm < RANSAC_distance;
distance_left_right = distance_left_right.*(double(myNorm));
[row, col] = find(distance_left_right);
ratio_threash=RANSAC_ratio;
sift_match_points_f_left=[];
sift_match_points_f_right=[];
for index = 1:numel(row)
    row_sort = sort(distance_left_right(row(index),:),'descend');
    ratio = row_sort(1)/row_sort(2);
    if ratio > ratio_threash
        sift_match_points_f_left = [sift_match_points_f_left ,f_left(1:2,row(index))];
        sift_match_points_f_right = [sift_match_points_f_right ,f_right(1:2,col(index))];
    end
end
% figure,showMatchedFeatures(img_left,img_right,sift_match_points_f_left',sift_match_points_f_right','montage');

%Step5 run RANSAC
% H = RANSAC(sift_match_points_f_left,sift_match_points_f_right);
max_homo=40;
hscores = zeros(1,max_homo);
% score homography
X1 = cat(1,sift_match_points_f_left,ones(1,size(sift_match_points_f_left,2)));
X2 = cat(1,sift_match_points_f_right,ones(1,size(sift_match_points_f_right,2)));
homographies = zeros(3,3,max_homo);
for ih = 1:max_homo
    points = randi(size(sift_match_points_f_left,2),1,4);
    A = [];
    b = [];
    for ip = 1:4
        point1 = sift_match_points_f_left(:,points(ip));
        point2 = sift_match_points_f_right(:,points(ip));
        p1_X = point1(1);
        p1_Y = point1(2);
        p2_X = point2(1);
        p2_Y = point2(2);
        A = cat(1, A, [-p1_X -p1_Y -1     0     0  0 p1_X*p2_X p1_Y*p2_X]);
        A = cat(1, A, [    0     0  0 -p1_X -p1_Y -1 p1_X*p2_Y p1_Y*p2_Y]);
        b = cat(1, b, [-p2_X; -p2_Y]);
    end    
    V = A\b;
%     [U,S,V] = svd(A);
    homographies(:,:,ih) = [V(1), V(2), V(3); V(4), V(5), V(6); V(7), V(8), 1];    
    X2_ = homographies(:,:,ih) * X1;
    du = X2_(1,:)./X2_(3,:) - X2(1,:)./X2(3,:) ;
    dv = X2_(2,:)./X2_(3,:) - X2(2,:)./X2(3,:) ;
    ok{ih} = (du.*du + dv.*dv) < 6*6 ;
    hscores(ih) = sum(ok{ih});
    X_X2_offset=X2_(1:2,:)-X2(1:2,:);
    avg_residuals(ih)=sqrt(sum(sum(X_X2_offset.*X_X2_offset)))/size(sift_match_points_f_left,2);
end
[value, index]=max(hscores);
residuals=avg_residuals(index);
inliers = value/size(sift_match_points_f_left,2);
judge_inlier=ok{index};
inliers_points_left=sift_match_points_f_left(:,judge_inlier);
inliers_points_right=sift_match_points_f_right(:,judge_inlier);

% figure,imshow(img_left,[0,255]);
% hold on;
% plot(inliers_points_left(1,:),inliers_points_left(2,:),'r*');
% figure,imshow(img_right,[0,255]);
% hold on;
% plot(inliers_points_right(1,:),inliers_points_right(2,:),'r*');

%Step6 wrap one image to the other
t = maketform('projective',homographies(:,:,index)');
% t = maketform('affine',homographies(:,:,index)');
[mosaic1 xdata1 ydata1] = imtransform(img_left_origin,t);
% imagesc(mosaic1), colormap gray;
% figure, imshow(mosaic1);

t_i = maketform('projective',eye(3));
% t_i = maketform('affine',eye(3));
[mosaic2 xdata2 ydata2] = imtransform(img_right_origin,t_i);
% imagesc(mosaic2), colormap gray;
% figure, imshow(mosaic2);

%Step7 create panorama
xMin = min(xdata1(1), xdata2(1));
xMax = max(xdata1(2), xdata2(2));
yMin = min(ydata1(1), ydata2(1));
yMax = max(ydata1(2), ydata2(2));

xRange = ceil(abs(xMin) + xMax);
yRange = ceil(abs(yMin) + yMax);
if size(size(img_left_origin),2)==3
    mosaic1_extend(:,:,1:3) = zeros(yRange,xRange,3);
    mosaic1_extend(1:size(mosaic1,1), 1:size(mosaic1,2),:)=mosaic1;
    mosaic2_extend(:,:,1:3) = zeros(yRange,xRange,3);
    mosaic2_extend(-round(yMin)+5:size(mosaic2,1)-round(yMin)-2, round(-xMin)+1:xRange-2,:) = mosaic2(2:size(mosaic2,1)-5,2:size(mosaic2,2)-1,:);
    mosaic=0.5*mosaic1_extend+0.5*mosaic2_extend;
    merge_img=mat2gray(mosaic);
    imshow(merge_img);
else
    mosaic1_extend(:,:) = zeros(yRange,xRange);
    mosaic1_extend(1:size(mosaic1,1), 1:size(mosaic1,2))=mosaic1;
    mosaic2_extend(:,:) = zeros(yRange,xRange);
    mosaic2_extend(-round(yMin)+5:size(mosaic2,1)-round(yMin)-2, round(-xMin)+1:xRange-2) = mosaic2(2:size(mosaic2,1)-5,2:size(mosaic2,2) - 1);
%     mosaic2_extend(-round(yMin)+5:size(mosaic2,1)-round(yMin)-2, round(-xMin)+1:xRange-2) = mosaic2(2:size(mosaic2,1)-5,2:size(mosaic2,2)-1);
    merge_img=0.5*mosaic1_extend+0.5*mosaic2_extend;
    imshow(merge_img);
end
end

